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Abstract 

We study the thermodynamic Casimir force between a spherical object and a plate. We consider 
the bulk universality class of the three-dimensional Ising model, which is relevant for experiments on 
binary mixtures. To this end, we simulate the improved Blume-Capel model. Following Hucht, we 
compute the force by integrating energy differences over the inverse temperature. We demonstrate 
that these energy differences can be computed efficiently by using a particular cluster algorithm. 
Our numerical results for strongly symmetry breaking boundary conditions are compared with the 
Derjaguin approximation for small distances and the small sphere expansion for large distances 
between the sphere and the plate. 
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I. INTRODUCTION 



In 1978 Fisher and de Gennes [Jj realized that when thermal fluctuations are restricted 
by the finite extension of the system, a force acts on the boundaries. Since this effect is 
similar to the Casimir effect, where the restriction of quantum fluctuations induces a force, 
it is called thermodynamic Casimir effect. Recently this force could be detected for vari- 
ous experimental systems and quantitative predictions could be obtained from Monte Carlo 
simulations of spin models 0]. The thermodynamic Casimir force is experimentally and 
maybe technologically interesting since it depends, in contrast to other forces, strongly on 
the temperature. Hence, since the temperature can be easily changed in experiments, the 
thermodynamic Casimir force can be easily manipulated. For a discussion of the thermody- 
namic Casimir force in the general context of nanoscale science see ref. 

Experiments on thin fllms of ^He and mixtures of ^He and "^He near the A-transition 
and on thin fllms of binary mixtures near the mixing-demixing transition 0, Isj have been 
performed. Theoretically these experiments are well described by a plate-plate geometry. 
For this geometry the thermodynamic Casimir force per area Fc follows the flnite size scaling 
form B 

Fc^kBTL-'d{t[L/^,,^Y/'') , (1) 

where L is the thickness of the fllm, t is the reduced temperature, ^o,+ is the amplitude 
of the correlation length in the high temperature phase, u is the critical exponent of the 
correlation length, d is the dimension of the bulk system and ^ is a universal function that 
depends on the bulk and boundary universality classes that characterize the fllm. Field 
theoretic methods do not allow to compute 9 for three-dimensional systems in the full range 



of the scaling variable [10|, yjj • Therefore the fact that 6 recently has been computed quite 
accuratel y by u sing Monte Carlo simulations of spin models constitutes valuable progress. 
In refs. (lUlS the three-dimensional XY bulk universality class and a vanishing fleld at 
the boundary have been studied, which is relevant for the experiments on ^He. A quite 
satisfactory agreement between the experimental results and the theory was found. In refs. 

the Ising bulk universality class and various types of boundary conditions were 



studied. Note that the mixing-demixing transition of binary mixtures belongs to the Ising 
bulk universality class. 

Experimentally it is easier to access the thermodynamic Casimir force for a sphere-sphere 
or a sphere-plate than for the plate-plate geometry. Several experiments for these geometries 
have been performed. In particular, the thermodynamic Casimir force between a colloidal 
)article immersed in a binary mixture of fluids and the substrate has been directly measured 



2J,|25|. In ref. |26| it was demonstrated that the thermodynamic Casimir forces in a colloidal 



system can be continuously tuned by the choice of the boundary conditions. In addition to 



homogeneous substrates also chemically patterned ones have been studied |27H29j . In other 



experiments, the thermodynamic Casimir force is the driving force for colloidal aggregation. 



see e.g. refs. |30|, |31 



For simplicity, we restrict the following discussion to the sphere-plate geometry. Further- 
more we assume that both the surfaces of the sphere and the plate are homogeneous. In 
this case, the thermodynamic Casimir force follows the scaling form j32| 

FciD,R,T) = ^K{t[D/i,^^Y'\A) (2) 

where R is the radius of the sphere, D the distance between the plate and the sphere and 
A = D/R. Note that in the literature often D/C,, where ^ is the bulk correlation length 
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in the high temperature phase, is used instead of t[D / ^q^^^/'^ as argument of the universal 
scahng function. The thermodynamic Casimir force is given by 



Fc{D^R^T) = -ksT^J^^ (3) 



where F{D,R, T) is the reduced free energy of the system. 

In the hmit D <^ R, the function Kjx, A) can be obtained from ^(t[L/^o,+]^'''^) by using 
the Derjaguin approximation (DA) l33|]. Instead, in the limit D ^ R the small sphere 
expansion |34] can be used. In ref. j32[ the function K is discussed for the case where both 
surfaces strongly adsorb the same component of the binary mixture. The behaviour of K 
for three-dimensional systems for D ~ i? is inferred from the exactly known one in two 
dimensions and from mean-field calculations. 

Here we discuss a method to determine the thermodynamic Casimir force for the sphere- 
plate geometry by using Monte Carlo simulations of spin models. The generalization to the 
sphere-sphere geometry is straightforward. The method is based on the geometric cluster 
algorithm of Heringa and Blote |35| . In the present work, we simulate the improved Blume- 
Capel on the simple cubic lattice. We focus on strongly symmetry breaking boundary 
conditions at the surfaces of the sphere and the plate. We assume that the surface of the 
plate prefers positive magnetisation. We study the two cases Sgphere = —1 and 1, where Sgphere 
is the value of the spins at the surface of the sphere. Preliminary results for Sgphere = 0, 
which are not discussed here in detail, show that also this case can be efficiently simulated. 
In our numerical tests we studied spheres of the radius R = 3.5,4.5,7.5 and 15.5 lattice 
spacings. For R = 3.5 and 4.5 we considered distances between the sphere and the plate 
up to D ~ lOR. We introduce an effective radius i?e// of the sphere which is determined 
by a particular finite size scaling method. To this end we compute the ratio of partition 
functions Z_ / by using a cluster algorithm which is closely related to the boundary fiip 
algorithm [36|, l37|. The index of Z gives the sign of s sphere- In this finite size scaling study 
we consider A f» 25, which is large compared with the values used otherwise. Analysing our 
data we find that replacing R by Rejf leads to a good scaling behaviour for the whole range 
of A studied here. As result we obtain reliable estimates for the scaling functions K^{x, A) 
and K+{x, A) for 1 ^ A ^ 12, where the subscript of K gives the sign of s sphere- The most 
striking observation is that in the low temperature phase /^-(x. A) deviates quite strongly 
from the Derjaguin approximation already for A ^ 1 and likely smaller. 

The present work should be seen as pilot study that allows for a number of improvements 
as discussed in section I VII The method discussed here should apply to a broader range of 
problems as we also shall discuss in section IVTl 

The paper is organised as follows. First we define the model. Then we discuss the 
geometry of the systems that we study. In section [IV] we discuss the algorithm that is used 
to compute the thermodynamic Casimir force. In section|V]we present our numerical results. 
First we discuss the performance of the algorithm. Then we give our results for the scaling 
functions K±{x,A). In section IVll we conclude and give an outlook. In the appendix we 
discuss the finite size scaling method that we used to determine the effective radius -Re//- 
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II. THE MODEL 



We simulate the Blume-Capel model on the simple cubic lattice. It is defined by the 
reduced Hamiltonian 

<xy> X X 

where G {—1, 0, 1} and x = (xq, Xi,X2) denotes a site on the simple cubic lattice, where Xi 
takes integer values and < > is a pair of nearest neighbours on the lattice. The partition 
function is given hj Z = ^^P(~-^)5 where the sum runs over all spin configurations. 

The parameter D controls the density of vacancies s^; = 0. In the limit D ^ —oo vacancies 
are completely suppressed and hence the spin- 1/2 Ising model is recovered. Here we consider 
a vanishing external field h = throughout. In d > 2 dimensions the model undergoes a 
continuous phase transition for — oo < D < D^^i at a /3c that depends on D. For D > Dt^i 
the model undergoes a first order phase transition, where Dfri = 2.0313(4) (38j . 

Numerically, using Monte Carlo simulations it has been shown that there is a point 
{D* , Pc{D*)) on the line of second order phase transitions^ where the amplitude of leading 
corrections to scaling vanishes. Our recent estimate is D* = 0.656(20) j39|. In 3^ we 



simulated the model at D = 0.655 close to Pc on lattices of a linear size up to L = 360. From 
a standard finite size scaling analysis of phenomenological couplings like the Binder cumulant 
we find /3c(0.655) = 0.387721735(25). Furthermore the amplitude of leading corrections to 
scaling is at least by a factor of 30 smaller than for the spin- 1/2 Ising model. In order to 
set the scale we use the second moment correlation length in the high temperature phase. 
Its amplitude is given by j4o| 



6nd,o,+ = 0.2283(1) - 1.8 X (z/ - 0.63002) + 275 x (/J^ - 0.387721735) 

using t = (3c — (3 as definition of the reduced temperature. (5) 

In the high temperature phase there is little difference between C,2nd and the exponential 
correlation length C. ^j, w hich is defined by the asymptotic decay of the two-point correlation 
function. Following 4ll |: 



lim^ = 1.000200(3) (6) 

for the thermodynamic limit of the three-dimensional system. Note that in the following 
always refers to ^2nd,o,+- 



III. GEOMETRY OF THE SYSTEM 

We have used fixed boundary conditions in 0-direction and periodic boundary conditions 
in 1 and 2 directions. In most of our simulations we have fixed s^, = 1 for xq = and 
Xq = Lq + 1. In order to check the effect of the second boundary, we have performed some 
simulations with = ^ for xq = and 3^ = for xq = Lq + 1. In both cases there are Lq 
layers with fluctuating spins. For i = 1,2 we take Xi = —Li/2 + 1,— Lj/2 -|- 2,..., Li/2 and 
Li = L2 = L throughout. 

To model the spherical object, we have fixed all spins Sx to Sgphere for 

^J{xo-hy + xj + xl < R (7) 
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Naively, the distance between the sphere and the plate is given by 



D = h-R . (8) 

Analysing our data we shall introduce an effective radius i?e// of the sphere that we deter- 
mine by using a finite size analysis in the appendix |X] below. Taking also into account the 
extrapolation length /g// at the boundary of the plate we arrive at 

Deff = h- R,ff + hff (9) 

where /g// = 0.96(2) — 0.5 = 0.46(2) for the model and boundary conditions discussed here 
(20| . Note that Ze// = 0.96(2) given in ref. 13] refers to Xq = 0.5 as location of the boundary. 



Mostly we have simulated the two = — 1 and 1. Here the fixed spins are coupled 

by the same /3 with their fluctuating neighbours as fluctuating spins and their fluctuating 
neighbours. One could imagine to chose the couplings depending on the position of the site 
on the surface of the sphere, in order to effectively improve the spherical shape of the object. 



IV. THE ALGORITHM 

In order to determine the thermodynamic Casimir force Fc{D, R, /3), eq. (E}, between the 
plate and the sphere, we compute finite differences of the reduced free energy F{D, R, (3) of 
the system 

- /3Fc(D, i?, /3) ^ AF{D, R, /3) = F{D + 1/2, R, /3) - F{D - 1/2, R, /3) . (10) 

Note that we chose the difference of the distances equal to one due to the discrete trans- 
lational invariance of the lattice. Similar to the proposal of Hucht [12^, we determine 
AF{D, R, (5) by integrating the corresponding difference of energies 

AF{D, R,(3) = - d(3 AE{D, R, (3) (11) 

where AE(D, R, (3) = E{D + 1/2, R, /5) - E{D - 1/2, R, (3) and 

E{D,R,^) = /j2'-'y) ■ (12) 

\<^y> I D,R,I3 

The value of (3o is chosen such that AE{D, R, /3o) ~ 0, which is the case for D ^ ibuikiPo)- 
The integration is done numerically, using the trapezoidal rule. It turns out that similar to 
the study of the plate-plate geometry O(IOO) nodes are needed to compute the thermody- 
namic Casimir force in the whole range of temperatures that is of interest to us. 

Naively one would carry out this programme by simulating the systems for D — 1/2 and 
D + 1/2 separately using standard Monte Carlo algorithms. In order to avoid systematic 
errors, Lq, L ^ D,R should be taken. However the variance of E{D, R, /3) and the CPU- 
time required for the simulations are increasing with the systems size. In fact, first numerical 
experiments have shown that even for moderate values of D and R, this approach fails to 
give sufficiently accurate results. 
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FIG. 1. Two-dimensional sketch of two spheres at xq = D — 1/2 and D + 1/2. The 0-axis is drawn 
in vertical direction. The overlapp of the two spheres is given by black squares with white dots. 
Sites that only belong to the lower sphere are represented by squares with vertical lines and those 
that only belong to the upper sphere by squares with horizontal lines. 



The expectation value < SxSy > only differs strongly between the two systems for nearest 
neighbour pairs < xy > close to the sphere and this difference decays at least power-like 
with increasing distance from the sphere. Therefore, computing AE{D, R, /3) one would like 
to concentrate on the part of the system, where the difference in < SxSy > is large and to 
avoid to pick up variance from the remainder. Naively one could implement this idea by 
ad hoc restricting the range of summation in eq. f|T2|) . The problem is that this ad hoc 
restriction leads to a systematic error that needs to be controlled, e.g. by comparing the 
results of different restrictions. In the following we shall discuss an algorithm, that generates 
the restriction of the summation to the neighbourhood of the sphere automatically and at 
the same time allows to compute AE{D, R, P) without bias. 



This algorithm is closely related to the geometric cluster algorithm [35|. Similar to ref. 



42[, two systems are updated jointly. These two systems share the values of /3,Lo,L and 
R. The distance of the sphere from the boundary is D — 1/2 and D + 1/2 for system 1 
and 2, respectively. In the following the spins are denoted by s^^i, where the second index 
indicates whether the spin belongs to system 1 or 2. Configurations are denoted by {s}i. 
The combined partition function is given by 

Zcom = 5^5^exp [-Hi{{s},) - H2{{S}2)] (13) 
where Hi is the reduced Hamiltonian of the system i. The cluster update comprises the 
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exchange of spins 



/ 



4,2 = Sx,l (14) 



for all sites x in the cluster. Performing this exchange of spins for all sites within a cluster is 
called flipping the cluster in the following. A cluster is build in the usual way: Links < xy > 



are deleted with the probability pd and frozen otherwise. In our case [35 



Pd = min [l,exp - 5^,2] Ki - 5^,2])] ■ (15) 

Two sites of the lattice belong to the same cluster if they are connected at least by one chain 
of frozen links. Various choices of the selection of the clusters to be updated are discussed in 
the literature. In the case of the Swendsen-Wang algorithm all clusters are constructed 
and then all spins within a given cluster are flipped with probability 1/2. In the single 
cluster algorithm ^] one first selects randomly one site of the lattice. Then one constructs 
only the cluster that contains this site. All spins within this cluster are flipped. Further 



choices have been discussed in the literature [45|, |46 . 

In our case, the choice of the clusters to be flipped is dictated by the requirement that 
the fixed spins have to keep their value under the update. The exchange of spins should 
not be allowed for those sites x, where Sx,i is fixed while Sx,2 is not, or vice versa. In the 
following these sites are called frozen sites. In Fig. [1] these sites are represented by squares 
with vertical or horizontal lines on them. In the following a cluster that contains at least one 
frozen sites is called frozen cluster. In the update all clusters except for the frozen ones are 
flipped. In order to perform this operation only the frozen clusters have to be constructed. 

Now let us discuss how AE, eq. (ITT]) , is computed. Let us assume that after equilibration 
we have generated a sequence of + 1 configurations that are labelled by t. We get 

1 ^ 1^1 
^ - ^ AE,st,t ~ ]y 5Z 2 + ^^-M+i) (16) 

t=i t=i 

where 

<xy> 

is the naive estimator of AE. 

Combining the right hand side of eq. f lT6|) and the particular properties of the update 
discussed above we arrive at the improved estimator 



AEimp - 2 {[Sx,2Sy,2 - Sx,lSy,l] + [4,24,2 ~ 4,l4,l]) 

<xy> 

= 2 {[Sx,2^y,2 - 4,l4,l] + [4,2-5^,2 " Sx,lSy^i]) 

<xy> 

= \ 5Z ([*^>2Sj;,2 - 4,lSy,l] + [s'x,2Sy,2 " Sx,lSy,l]) (18) 

<xy>&Cf 

where < xy >G Cf means that at least one of the sites x or y belongs to a frozen cluster. 
Note that for our choice of the update 

sL,i4,i = ^^,2Sy,2 and s^_24,2 = Sx,iSy,i (19) 
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for all nearest neighbour pairs < x,y > where neither x nor y belongs to a frozen cluster. 

Whether we have achieved progress this way depends on the average volume taken by 
frozen clusters. Our numerical study discussed in detail below shows indeed that for the 
whole range of /3 the total size of all frozen clusters is small compared with the volume of 
the lattice. Hence the variance of AEi^p is drastically reduced compared with the naive 
estimator. Furthermore, since the sum runs only over the frozen clusters, the numerical 
effort to compute AEi^p is much smaller than the one to compute the naive estimator. 

In the following we shall discuss the details of our implementation. We simulate a large 
number of distances Dmin, Dmin + ^,---,Dmax jointly. We perform the exchange cluster 
update and the measurement of AEimp for all Dmax — Dmin pairs {Di, D2) with D2 — Di = 
1. In our C-program we store the spins in the array int spins [ID] [LO] [L] [L] ; where 
ID= Dmax — Dmin + 1- In Order to save CPU time we do not copy all spins outside the 
frozen clusters from spins [ic] [] [] [] to spins [ic+1] [] [] [] and vice versa. Instead, we 
do that for the spins that belong to frozen clusters. This way, the systems with the distances 
Di and D2 of the sphere from the plate interchange their position in the array spins. In 
order to keep track of where the systems are stored in the array spins, we introduce the 
array int posi[ID] ; where the index is given by id= D — Dmin and posi[id] indicates 
where the system with the distance D of the sphere from the plate is stored in spins. 
Note that parallel tempering simulations are usually organized in a similar way. Instead of 
swapping configurations, the location of the systems with the temperatures Tj and Tj+i is 
interchanged. 

In order to get an ergodic algorithm we supplement the exchange cluster algorithm dis- 
cussed above with standard updates of the individual systems. To this end we perform 
Swendsen-Wang [43[ updates and heat-bath updates. In case of the Swendsen-Wang cluster- 
update we keep the clusters that contain fixed spins fixed. Those clusters that do not contain 
frozen spins are flipped with probability 1/2. 

These updates are performed in a fixed sequence: First we sweep once for each systems 
over the whole lattice using the heat-bath algorithm. It follows one Swendsen-Wang cluster- 
update for each system. Then we performed ^m-times the following two steps: 

• For D = Dmin, up to Dmax — 1 perform in a sequence the exchange cluster update 
and the measurement of AEimp for the pairs of distances D and D + 1. 

• Perform for each system one sweep with the heat-bath algorithm over a sub-lattice 
characterized by xq < Iq, — ^i/2 < xi < li/2 and — ^2/2 < X2 < ^2/2. 

In the following we shall denote this sequence of updates as one cycle of the algorithm. Since 
the estimator AE^mp takes mainly data from the neighbourhood of the sphere, it seems useful 
to perform additional updates in this region of the lattice. Therefore we introduced the heat- 
bath updates of the sub-lattice. Here we made no attempt to optimize the shape and size 
of this sub-lattice nor the frequency of these updates. 

V. NUMERICAL RESULTS 

As a first test of our implementation we simulated systems with the radius R = 2.5 and 
lattices of the size 28 x 16^ and 38 x 20^. We considered the distances 1.5 < D < 9.5 between 
the sphere and the plate. We fixed Sx = I for xq = and xq = Lo + 1. We simulated both 
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Ssphere = "1 and +1 at the two values 13 = 0.387 and 0.389 of the inverse temperature. 
We simulated with the exchange cluster algorithm and without. For these lattice sizes 
we still get differences of the energies AE that are clearly larger than the statistical error 
from simulation that take about one day of CPU-time on a single core of the CPU in the 
naive way. Our results with and without the exchange cluster algorithm are compatible 
within the statistical errors. But even for the small lattices sizes used here, the statistical 
error, using an equal amount of CPU-time, is reduced by more than a factor of 10 by using 
the exchange cluster algorithm and the improved estimator f|T8|) . After these preliminary 
studies were concluded successfully, we performed an extensive study aiming at physically 
relevant results. First we summarize the parameters of our simulations. Then we discuss 
the performance of the algorithm, focussing on the size Ndu of the frozen clusters minus the 
number of fixed spins. Finally we present our results for the thermodynamic Casimir force 
and compare them with the literature. 

A. Our simulations 

We studied spheres of the radii R = 3.5, 4.5 and 7.5 for both Sgphere = —1 and 1 and 
R = 15.5 for Ssphere = 1 ouly. For these radii, we simulated systems with 1.5 < D < Dmax, 
where D^nax = 16.5, 40.5, 32.5, and 10.5 for R = 3.5, 4.5, 7.5 and 15.5, respectively. In the 
following we shall denote the number of values of D by Ijj. We simulated these systems 
for about 100 values of (3. We adapted the step-size in /3. Close to f3c, the step-size is the 
smallest. For all values of /3 we simulated lattices of the size 148 x 80^, 298 x 100^, 298 x 160^ 
and 398 x 200^, respectively. In the case of these simulations, we fixed Sx = I for xq = 
and xo = Lq + 1. We fixed the parameter of the update cycle to im = 20. The sizes of 
the sub-lattice that is updated in excess im times per cycle using the heat-bath algorithm is 
28 X 16^, 55 X 20^, 56 x 32^ and 58 x 64^ for R = 3.5, 4.5, 7.5 and 15.5, respectively. These 
parameters are chosen ad hoc. For lack of human time we did not try to optimize them. 

In order to check for the effect of the finite values of Lq and L, we redid the simulations 
for R = 3.5 in the neighbourhood of (3c on lattices of the sizes 198 x 120^ and 298 x 200^, 
where we considered both = and = 1 for xq = Lq + 1. For R = 4.5 we simulated 
in addition lattices of the size 398 x 150^ with 3^ = for xq = Lq + 1 and 448 x 150^ with 
Sa, = 1 for Xq = Lq + I in the neighbourhood of Pc- For these additional simulations we chose 
im = 10. The sizes of the sub-lattice that is updated in excess are the same as above. In 
order to keep the figures readable we shall give in sections IV C\ IV Dl and IV El below results 
for s^^ = 1 at Xo = i^o + 1 only. The results obtained for s^^ = at Xq = -^^o + 1 confirm our 
final conclusions. 

Focussing on larger distances we simulated R = 3.5 for 25.5 < D < 32.5. For Sgphere = —1 
in the range 0.3864 < (3 < 0.3888 for Ssphere = 1 in the range 0.3864 < (3 < 0.3892 we 
simulated a lattice of the size 598 x 200^. The size of the sub- lattice that is updated in 
excess im times per cycle using the heat-bath algorithm is 43 x 16^. For (3 less close to (3c 
we simulated smaller lattices. For these simulations we took im = '20. Finally we performed 
simulations for R = 4.5 at the distances 31.5 < D < 40.5 using a lattice of the size 598 x 200^ 
in the range 0.386 < (3 < 0.3886 also taking im = 20. 

For each of these simulations we performed 10^ up to 2 x 10^ update cycles. In all 
our simulations we used the Mersenne twister algorithm ji^] as pseudo-random number 
generator. Our code is written in C and we used the Intel compiler. Simulating at (3 = (3c, 
our Swenden-Wang update takes about 1.1 x 10~^ s and the local heat-bath update about 
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FIG. 2. We plot the average size Ndu of the frozen clusters minus the number of frozen sites for 
the radius R = 7.5 and D = 32 as a function of the inverse temperature /3. We give the results for 
s sphere = " 1 and 1. The dotted line should only guide the eye. 



3.5 X 10^* s per site on a single core of a Quad-Core AMD Opteron(tm) 2378 CPU. The 
exchange cluster, including the measurement ( JT8l) takes roughly 2.6 x 10^'' s per site of the 
frozen clusters. In total our simulations took roughly the equivalent of 50 years of CPU time 
on one core of a Quad-Core AMD Opteron(tm) 2378 CPU. 



B. Size of the frozen clusters 

The CPU-time needed for one update cycle depends on the size of the frozen clusters. 
Also the variance of the estimator ( fTSi) compared with the standard one can be small only 
if the size of the frozen clusters is small compared with the volume of the lattice. Therefore 
we shall discuss the size of the frozen clusters and its dependence on the parameters /3, R, 
D and the lattice size in some detail below. 

To give a first impression, we plot in Fig. [2] the average size N^iu of the frozen clusters 
minus the number of frozen sites for the radius R = 7.5 and the distance D = 32, i.e. for 
the pair of systems with D = 31.5 and 32.5, for Sgphere = —1 and 1. For decreasing {3 in 
the high temperature phase, the average size of the frozen clusters for s sphere = — 1 and 1 
seems to become identical. As the critical point is approached, the size for s sphere = — 1 
becomes clearly larger than for s sphere = 1- This remains true in the low temperature phase. 
In both cases, the maximum Nciu,max is located close to (3c- We find Nciu,max ~ 3900 and 
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FIG. 3. We plot Nciu as a function of /3 for i? = 3.5, D = IQ and Sgphere = 1- We focus on the 
neighbourhood of the critical point. We give results for the lattices sizes 148 x 80^, 198 x 120^ and 
298 X 200^. 



Nciu,max ~ 1190 for Sgphere = ~1 and 1, respectively. Note that this is only a small fraction 
of the total volume of the lattice, 298 x 160^ = 7628800. 

Let us now discuss the behaviour of Ndu in more detail. Since = 1 for /3 = 0, the 
frozen clusters consist of the frozen sites only, i.e. Ndu = 0. At low temperatures, the 
configurations become ordered. Therefore also in this limit the frozen clusters consist of the 
frozen sites only. 

From our numerical data we see that for fixed /3 and lattice size, Ndu is monotonically 
increasing with increasing distance D. For Lq, L,D ^ C,, Ndu seems to converge to a finite 
value Ndu,oo- In the high temperature phase, (3 < (3c, this value is the same for Sgphere = — 1 
and 1. We fitted our data in the high temperature phase with the ansatz Ndu,oo = c{(3c—(3)~^ . 
For R = 3.5 we get quite reasonable fits with the result y ~ 0.6. Fits for larger radii, where 
we have less data apparently give smaller values of y. However, since Ndu,oo should increase 
with increasing radius i?, this observation seems to be an artifact of our finite data set. 
Here a theoretical understanding of the behaviour of Ndu,oo, relating y with known critical 
exponents would be desirable. 

Next let us study the dependence of Ndu on the linear lattice sizes Lq and L in the 
neighbourhood of the critical point. In Fig. Owe plot Ndu as a function of (3 for R = 3.5, 
D = 16 and Ssphere = 1 for the lattice sizes 148 x 80^, 198 x 120^ and 298 x 200^. For 
(3 ^ 0.387 and /3 ^ 0.3882 the estimates of Ndu for the different lattice sizes are consistent 
within the statistical error. Instead for 0.387 ~ /3 ~ 0.3882 a dependence of Ndu on the 
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FIG. 4. We plot Nciu as a function of D ior R = 4.5, Ssphere = —1 and /3 = /3c. We give results for 
the lattice sizes 298 x 100^, 448 x 150^ and 598 x 2002. 



TABLE I. 7Vc/« for 

s sphere — ~1 3,iid 1 at the Critical point for D A.3R. The results are taken 
from the lattice sizes 148 x 80^, 298 x 100^, 298 x 160^ for R = 3.5, 4.5 and 7.5, respectively. 



R 


^sphere 1 


^sphere 1 


3.5 


1096(3) 


269(1) 


4.5 


2464(6) 


565(2) 


7.5 


3862(8) 


970(2) 



lattice size can be seen. For most values of (3 in this range N^u is increasing with increasing 
lattice size. However for /3 ^ 0.388 the value of N^u for 148 x 80^ seems to be a bit larger 
than those for the lattices sizes 198 x 120^ and 298 x 200^. For Sgphere = ~1 we get a 
qualitatively similar behaviour. The same holds for R = 4.5 for both s sphere = ~1 and 1. 

Finally let us discuss the behaviour of Ndu at the critical point in more detail. In Fig. 
IHwe plot Nciu as a function of D ioi R = 4.5, s sphere = — 1 and (3 = f3c. The value of Ndu 
seems to increase somewhat faster than linearly with the distance D. Differences between 
the results for the lattice sizes 298 x 100^, 448 x 150^ and 598 x 200^ are clearly visible. Similar 
observations can be made for Sgphere = 1 and also for other the radii that we simulated. 

In table [I we give N^iu at (3c for the radii R = 3.5, 4.5 and 7.5 for the distance D ^ 4.3i?, 
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i.e. D = 15, 19, and 32, respectively. Comparing the results for R = 4.5 and 7.5 one 
might conclude that Ndu is increasing less than linearly with the radius R, which is quite 
encouraging with respect to the scaling of the performance of the algorithm. 

In section [IV] we argued that the estimator (fTSjl reduces the variance if Ndu ^ Lq x L^. 
This is indeed the case for the full range of /3 and all radii and lattice sizes studied here. A 
better theoretical understanding of the behaviour of Ndu is of course desirable. 



C. The thermodynamic Casimir force: Effects of the finite lattice 

For P ^ Pc corrections to the limit Lq, L — )■ oo decay exponentially with increasing lattice 
size. Therefore corrections should be negligible for Lq, L ^ ibuik- Instead, at the critical 
point, we expect that corrections decay power-like. For R = 3.5 we checked explicitly how 
AE depends on the lattice size. To this end we simulated in the neighbourhood of Pc in 
addition to the lattice of the size 148 x 80^ ones of the sizes 198 x 120^ and 298 x 200^. 
In Fig. [5] we plot for D = IQ our results for the difference of energies AE as a function 
of p. We give our results for s sphere = — 1 and 1 in the upper and the low part of Fig. |5l 
respectively. We see differences between the results obtained for the 148 x 80^ and the larger 
lattices that are statistically significant for 0.3868 ~ /3 ~ 0.3886 and 0.3868 ~ /? ~ 0.3882 
for Ssphere = "1 and 1, respectively. Note that ^b„/fc(0.3868) ~ 18.66, using eq. 

In Fig. [6]we give the numerical results of —AF for R = 3.5, D = IQ and Sgphere = —1- In 
the upper part of the figure we give the result of the numerical integration of AE obtained 
by using the data for the 148 x 80^ lattice. In the lower part of the figure we focus on the 
neighbourhood of Pc, where effects due to finite Lq and L become important. As check, we 
computed — AF in an alternative way: We replaced the estimates of AE from the 148 x 80^ 
lattice, where available, by those from the 198 x 120^ or the 298 x 200^ lattice. We find that 
the behaviour of —AF is changed only little. The difference between these results is of a 
similar size as the statistical error. Note that the statistical error of — AF is the accumulated 
error of AE over the range of the integration. 

We conclude that in a neighbourhood of the critical point the effect of the finite lattice 
size on AE is clearly visible. However, since the range of P where this is the case, is quite 
small, — AF is affected less. The amplitude of this effect can be estimated by comparing 
the results obtained for different lattice sizes. 



D. The scahng function K^{x,A) for Sgphere = 1 

In this section we discuss our results for the scaling function A) for s sphere = 1- In a 

first step we check how well our results obtained for different radii collapse on a single curve. 
To this end we make use of the effective radius determined in appendix \^ see table El and 
the effective distance ([9]). In particular A = Deff/Reff in the following analysis of our data. 
In the plots given below we show the statistical error of the data only. The uncertainty of 
Deff and -Re// still comes on top of this. Essentially it is a two percent uncertainty of the 
scales. 

In figure [7] we plot -R^ffAFA'^ function of x = tlD^ffm^l" for A ^ 2.7 and 
4.8. We multiplied by A^ to reduce the dependence on A. Note that for small values of 
A, /s:(0. A) oc A-2 and large ones /s:(0. A) oc A-^'"-^ IH, Q, where P/v = (1 + r/)/2 



and T] = 0.03627(10) We took the data obtained by using 148 x 80^, 298 x 100^ and 
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FIG. 5. Numerical results Ai? for the radius R = 3.5 and the distance D = 16. The lattice sizes 
are 148 x 80^ (black circle), 198 x 120^ (blue upward triangle), and 298 x 200^ (red downward 
triangle). In the upper part of the figure we give results for s sphere = ~1 and in the lower part 

results for S sphere = 1- 
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FIG. 6. Numerical results for —AF for the radius R = 3.5, the distance D = 16 and Sgphere = — 1- 

In the upper part of the figure we give the result obtained by using a 148 X 80^ lattice for the full 

range of /3 that we considered. In the lower part of the figure we focus on the neighbourhood of /3c- 

We plot the results for the lattice sizes 148 x 80^ (black circle), 198 x 120^ (blue upward triangle), 

and 298 x 200^ (red downward triangle). For a discussion see the text. 
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298 X 160^ lattices for the radii R = 3.5, 4.5 and 7.5, respectively. For this choice L/R is 
roughly the same for the three radii considered. 

In particular for A 4.8 we see an excellent collapse of the data obtained for the three 
different radii. Note that things look far worse when R and D are used instead of -Re// and 
Deff. For A ^ 2.7 the collapse of the data is slightly worse than for A ^ 4.8. In particular 
for X ^ 10 we see some deviation of the result for R = 3.5 from that for the two other radii. 
This is not too surprising, since A ^ 2.7 means D = 6 for R = 3.5 and also the values of f3 
that correspond to a; ^ 10 are far from Pc- 

We conclude that the data scale well and our numerical results should describe the scaling 
limit. At the level of our accuracy, one should interpret data obtained for ^ 6 with caution. 

In Fig. IS] we plot i^+(a;, A)A^ as a function of a; = t[D/^QY^'^ for various values of A in 
the range 0.75 < A < 10.89. For each value of A, we took the largest radius, where data 
are available. For comparison we plotted the Derjaguin approximation in the limit A — )■ 0. 



For a brief discussion of the Derjaguin approximation see e.g. ref. [32|. As input we used 
the results for 6*++ (a;) and 6+-{x) obtained in ref. (23| . 

We find that the shape of Kj^{x, A)A^ changes only little with varying A. In particular 
the location of the minimum changes very little with varying A. The absolute value of 
K^{x, A)A^ is increasin g wi th increasing A, which is consistent with the fact that for large 
A, K+{0,A) oc A-^/'^-i 

In the following we study in more detail how the minimum of K^{x,A) as function of 
X depends of A. Fortunately, Xmin is large enough, such that finite Lq and L effects are 
well under control. We located the minimum of —AF by finding the zero of AE. To 
this end, we interpolated AE linearly in (3. Using these results, we computed x^m = 
iPc — l3min)[Def f / ^o]^^'^ ■ In Fig. [9] we plot Xmin as a function of A using our data obtained 
for the radii R = 3.5, 4.5, 7.5 and 15.5. We see that for A ^ 2 the results obtained for 
different radii agree quite nicely. At A = 2 the estimate Xmin ~ 2.43 seems plausible. With 
further increasing A the value of Xmin is slightly decreasing. For A ^ 12 we get Xmin ~ 2.3. 
For A ^ 2 clear deviations between the results obtained by using different radii R can be 
observed. Therefore we abstain to give a conclusion for the scaling limit in this range. Note 
that the Derjaguin approximation gives Xmin ~ 2.52 for the limit A — )■ 0. The Derjaguin 
approximation predicts that Xmin is increasing with increasing distance A. If that prediction 
is correct, Xmin should show a maximum in the interval < A ^ 1. 

Our results can be compared with the upper part of FIG. 2 of ref. j32|. Note that the 
authors of ref. j32| consider D/C, ^ x'^ as variable. They obtain Xmin = l.Q^/^-^^^^"^ = 2.1... for 
A = using the Derjaguin approximation. This discrepancy should be due to the different 
numerical estimates of the scaling function 6'++ of the plate-plate geometry. E.g. following 
FIG 3 of ref. (slj ^++(0) = -0.326, while our present estimate is ^++(0) = -0.410(8) |23|. 
Then, according to ref. j32| and reaches a maximum at A ^ 3 assuming 

the value Xmin ~ 2.5^/°-^3°°2 ^ ^ ^i^g^^^ ^j^^ increasing A it is decreasing and reaches 
Xmin ~ 1.51/0-63002 ^ -^ g limit A — )■ oo. In particular for the value assumed at A ?a 3 

we see a clear discrepancy with our results. It is beyond the scope of the present work to 
discuss the reliability of the first order e-expansion that has been used in ref. |32| to get the 
small-sphere expansion. 

In Fig. [TU] we plot i^+(xmm,^)^^ as a function of A. We find a nice collapse of the 
data obtained from different radii of the sphere. In contrast to Xmin, this is also the case 
for A ^ 2. Taking into account our Monte Carlo data only, K+{xmin, A)A'^ seems to be a 
monotonically decreasing function of A. However we note that for A ^ 0.4 our numerical 
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FIG. 7. For 

^sphere — 1 ^6 plot — R^f f /^F /S."^ as a function of x — t[D^j j- /^^q]^^'^ the radii R — 3.5, 
4.5 and R = 7.5. In the upper part we plot results for A ~ 2.7, which means D = 6, 9 and 18 
for R = 3.5, 4.5 and 7.5, respectively. In the upper part we plot results for A ~ 4.8, which means 
D = 12, 17 and 32 for R = 3.5, 4.5 and 7.5, respectively. Note that A = D^ff/Reff is used. 
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FIG. 8. We compare the Derjaguin result for lirriA-s-o ^)^^ with our Monte Carlo data for 
-fC+(x, A)A^ for various values of A. In particular we have used our data for K = 15.5 and D = 10 
giving A ^ 0.75, R = 7.5 and D = 11, 18, 25, and 32 giving A ^ 1.76, 2.78, 3.80, and 4.83, 
respectively, and R = 4.5 and D = 21, 29 and 40 giving A 5.87, 7.98, and 10.89, respectively. 



estimates of K^{xmin, A)A^ are overshooting limA->.o K+{xmin, A)A^ ^ —4.63 obtained from 
the Derjaguin approximation. Also note that the Derjaguin approximation predicts that 
i^+(a;mm, A)A^ increases with increasing A, which is indicated by a dashed black line in 
Fig. IHl Therefore i^+(xmm, A)A^ should show a maximum in the interval < A ^ 0.4. In 
the lower part of Fig. 2. of ref. j32[ /('^(xmmj A)A^+^/'^ is shown as a function of A. For 
A 10 it shows a maximum (minimum with the sign of ref. j32| ) taking a value slightly 
larger than —3. Taking our Monte Carlo data, K^{xmin, A)A^"'"^/'^ is still slightly increasing 
at A ^ 10 and takes the value ~ —3.25. 

Finally we discuss the scaling function at the critical point, i.e. x = 0. In figure [TT] we 
plot our numerical results for -ft'+(0, A)A^ as a function of A. At the critical point we have 
to consider finite L and Lq effects. Therefore, for R = 4.5 we plotted our results obtained 
by using lattices of the sizes 298 x 100^, 448 x 150^ and 598 x 200^. In particular for large 
A we see a clear deviation of the results obtained for the 298 x 100^ lattice from the other 
two. On the other hand the estimates obtained for the 448 x 150^ and 598 x 200^ lattices 
are consistent within the statistical error. Therefore deviations from the Lo,L — )■ oo limit 
should be of similar size as the statistical error. In the case of R = 3.5 we plot the data 
for the largest lattice size available. We find a quite good collapse of the data obtained for 
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FIG. 9. We plot our numerical estimates of Xmin as a function of A 



different radii. For A = 1 we get K+{0,A)A^ ^ -2.8. We see that K+{0,A)A^ is almost 
linearly decreasing with increasing A. For A = 11 we read off K^{0, A)A'^ ^ —6. The 
last three data points for R = 3.5 seem to suggest that -ft'+(0, A)A^ is increasing again for 
A ^ 11. However this is not statistically significant and is likely an artifact. 

Using the Derjaguin approximation, we get limA-s-o -^+(0, A)A^ ^ —2.59 and as for 
K+^min, -^+(0,A)A^ is increasing with increasing A. Therefore it seems plausible that 
-ft'+(0, A) A^ has maximum in the range < A ^ 1. 



E. The scaling function K-{x,A) for Sgphere = —1 

Next let us turn the scaling function K_{x,A) for s sphere = —1- In figure [12] we plot 
i^_(x,A)A^ as a function of x for A — )■ and A 1. The limit A — )■ is obtained 
by using the Derjaguin approximation. Unfortunately we do not have data for the scaling 
function for the plate-plate geometry for x ^ —48 and the amplitude of 6'+_(x) in 

this range is still significant. Therefore we have extrapolated in two different ways: 

(1) 9^ (x) = const for x ^ —48 and (2) 9^ (x) is linear in the interval —65 < x < —48 

and e+_(-48) = 0.8, 9+-{-Qb) = 0, for x < -65 we take 9+-{x) = 0. The true 9+-{x) 
should be enclosed by these two extrapolations. In figure [121 the results based on these two 
extrapolations are denoted by Derjaguin 1 and Derjaguin 2, respectively. For comparison 
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FIG. 10. We plot our numerical estimates of A)A^ as a function of A. 



we plot our Monte Carlo result obtained for R = 7.5 and D = Q which corresponds to 
A ^ 1. We observe that the two curves are very similar in the high temperature phase. 

Since 6^ {x) has its maximum in the low temperature phase, the same holds for K_{x,A) 

computed by using the Derjaguin approximation. In contrast, for A «i 1 we find that the 
maximum is located in the high temperature phase, even though very close to the critical 
point. Furthermore, K_{x, A)A'^ at A = 1 is much smaller than lim^^Q K^{x, A)A'^ in the 
low temperature phase. Also the decay of K^{x, A)A^ as x — )■ — oo seems to be much faster 
for A f» 1 than for A — )■ 0. 

This observation might be explained as follows: Deep in the low temperature phase the 
physics of a film with H — boundary conditions is governed by the interface between the 
two phases of positive and negative magnetisation. In the case of the film, this interface 
can move quite freely between the two plates. In contrast, in the case of the sphere-plate 
geometry, the interface is closely attached to the sphere, minimizing the area of the interface. 

In figure [T^ we plot K_{x, A)A^ as a function of x for various values of A. With increasing 
A the maximum of K-{x, A)A^ increases and also Xmax slowly increases. For the range of A 
studied here, K^{x, A)A^ is increasing with increasing A for x ^ —0.6, while it is decreasing 
with increasing A for x ^ —0.6. 

Next we study in more detail the maximum of K-{x,A) as a function of x. In figure 
[14] we plot the location Xmax of the maximum as a function of A. Since Xmax is very close 
to zero, finite Lq and L effects might be large. In order to check for these effects, we plot 
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FIG. 11. We plot our numerical estimates of i^-(-(0, A)A^ as a function of A. 



for R = 4.5 the results obtained by using lattices of the sizes 298 x 100^, 448 x 150^ and 
598 X 200^. In fact for A ^ 3 we see that the results obtained for the 298 x 100^ lattice 
clearly deviate from those for the other two lattice sizes. On the other hand, the results for 
the lattice sizes 448 x 150^ and 598 x 200^ are essentially consistent and are more or less 
consistent with those for R = 3.5. Therefore we consider the results that we have obtained 
for our largest lattice sizes as reasonable approximation of the Lq, L — )■ oo limit. 

Using the Derjaguin approximation, we get Xmax ~ —1.43 in the limit A — )■ 0. Within 
the Derjaguin approximation Xmax is decreasing with increasing A. In contrast for all values 
of A that we studied by Monte Carlo simulations, we find Xmax > 0. Looking at our data 
obtained for R = 7.5 we estimate x^az = at A fti 0.5. x^ax slowly increases with increasing 
A. For A ^ 12 we get Xmax ~ 0.6. 

In figure [T5] we plot i^_(xmax, A)A^ as a function of A. Also here we check for finite 
Lq and L effects by plotting for _R = 4.5 the results obtained for three different lattice 
sizes. Also here we find that the result for the 298 x 100^ lattice differs from those for the 
448 X 150^ and 598 x 200^ lattices for larger values of A, while the results obtained for 
lattices of the sizes 448 x 150^ and 598 x 200^ are in quite good agreement. For A ^ 4 the 
results for R = 7.5 and the other radii do not agree within error bars indicating violations of 
scaling that are larger than the numerical errors. Taking into account these observations we 
conclude that i^'_(xmax, A)A^ increases roughly linearly in the range of A that we studied. 
We get K_{x„^ax, AjA^ ^ 16 for A ^ 1 and K_{xmax, A)A2 ^ 21.5 for A ^ 12. In the limit 
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FIG. 12. We plot K_{x, AjA^ as a function of x for A ^ and A 1. 



A — >• we get K_{xmax, A)A^ ^ 19.2. Within the Derjaguin approximation K_{xmax, A)A^ 
is decreasing with increasing A. Hence i^'_(xmax, A)A^ should reach a minimum in the 
interval < A ^ 1. 

In figure [TBI we plot A'_(0, A)A^ as a function of A. Our observations concerning finite 
Lq and L effects and scaling are qualitatively similar as for K_(xmaa:, A)A^. However the 
amplitudes of the differences are larger than for i^_(xmaa;, A)A^. For large values of A 
even the results for the lattice sizes 448 x 150^ and 598 x 200^ do not agree within the 
statistical error. Taking into account these observation we conclude that -ft'_(0,A)A^ is 
increasing roughly linearly in the range 1 ^ A ^ 12. We get -ft'_(xmax; A)A^ ^ 16 for 
A fa 1 and -ft'_(xmax; A)A^ ^ 19.5 for A ^ 12. In the Derjaguin approximation we get 
limA-s.0 -^-(0, A)A^ ^ 17.6. In the Derjaguin approximation, i^_(0,A)A^ is decreasing 
with increasing A. Therefore, as it is the case for K_{xmax, A)A^, K_{0, A)A^ should show 
a minimum in the interval < A ^ 1. 



F. Small sphere expansion at the critical point 



In the limit D ^ R the behaviour of the reduced excess free energy at the critical point 
is given by |3J] 



FUD,R) = -^ —] (20) 
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where is the dimension of the field ip. The amphtudes are defined by 

{i^oA)buik = B^r-^-^ (21) 

and 

{i^r)Usr>ace = At{2r)-^^ (22) 

where in the case of eq. f l2^ r is the distance from the boundary of type a. For the boundary 
conditions discussed here, the leading contribution comes from the order parameter with 

x^ = /3/u=^{d-2 + 7]) (23) 

where 7] = 0.03627(10) j39|. The subleading contribution is due to the energy density e with 

x, = d-l/u (24) 

where u = 0.63002(10) H. 

In relation with ref. ^J0\ we had simulated 512'^ lattices with ++ and +0 boundary condi- 
tions in one direction and periodic ones in the other two directions at /3 = 0.387721735(25). 
We had measured the magnetisation profile. Taking into account the extrapolation length, 
we get from the magnetisation profile at r ^ 30 the estimate 

AI = 1.1268(5) . (25) 



23 



sphere 




FIG. 14. We plot 

Xmax 8'S 3- function of A for Sgph^^e — ~1- We give our results for the radii 
R = 3.5, 4.5 and 7.5. In the case of R = 4.5 we compare the results obtained for three different 
lattice sizes. 



In order to compute the amplitude we simulated systems with periodic boundary con- 
ditions in all three directions. We simulated lattices of the sizes 128^, 256'^, 512^ and 1024^ 
at 13 = 0.387721735(25). We measured {s^Sy) for y - x = (i,0,0), (0,i,0) and (0,0, and 
i = 1, 2, ... . Using the numerical results of the two-point function we computed 

B<t>,eff{l) = ^^''^(SxSx+(i,0,0)) (26) 

using 2x^ = 1.03627(10). It turns out that B^^eff{i) shows a shallow minimum for a rather 
small value oii. At larger distances, B^^(.ff{i) increases due to the periodic boundary condi- 
tions imposed. The minimum is given by {imin, B^^eff,min) = (4, 0.19397(8)), (6, 0.19026(15)), 
(7,0.18838(11)), and (8,0.18745(13)) for lattices of the size 128^, 256^, 512^ and 1024^, re- 
spectively. Assuming power like corrections, we arrive at our final estimate 

5,^ = 0.1865(3) . (27) 

Combining eqs. fl25|27p we arrive at 

a = ll±il± = 6.808(12) (28) 

B^ 
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FIG. 15. We plot K^{xmax, ^)^'^ as a function of A for Sgphere = —1- We give our results for 
the radii R = 3.5, 4.5 and 7.5. In the case oi R = 4.5 we compare the results obtained for three 
different lattice sizes. 



This result can be compared with a ~ 7.73 given in ref. 32|, where the authors had 



deduced the estimate for three dimensions from the exact result in two dimensions and the 



e-expansion. For a discussion see footnote [13] of ref. [32 . 

As check we computed for the radii i? = 3.5 using the sphere flip algorithm 

discussed in the Appendix. We have simulated lattices of the size Lq x where we have 
fixed = 1 for xq = and s^, = for xq = Lq + 1. Our numerical results for R = 3.5 are 
summarized in table [TTl First we have simulated for h = D + R = 72 lattices of different 
linear sizes Lq and L in order to check numerically the effect of the finite size of the lattice. 
We expect that at the critical point, the result for converges with a power law to the 

Lo, L cx) limit. The resuhs of (Lq, L) = (499, 400) and (Lq, L) = (799, 600) are consistent 
within the statistical error. Next we simulated lattices of the size {Lo,L) = (499,400) for 
various values of h = D + R. 

Inserting = —At into eq. fl20l) we get ioi D ^ R 



- ln(Z_/Z,) = 2^ ( A . (29) 
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FIG. 16. We plot K^{0, A)A^ as a function of A for Sgphere = —1- We give our results for the radii 
R = 3.5, 4.5 and 7.5. In the case of R = 4.5 we compare the results obtained for three different 
lattice sizes. 



For finite D we define 

where we replaced R and D by Rejf, table El and D^fj, eq. (jHl), respectively, to reduce 
corrections to scaling. 

In the last column of table [III we give estimates of a obtained for various distances D. It 
is decreasing with increasing distances. For the largest distance that we simulated, it is still 
clearly larger than the estimate (l28l) . indicating that still at A 30 subleading contributions 
to Fex are significant. 



VI. CONCLUSIONS AND OUTLOOK 



We have demonstrated how the thermodynamic Casimir force for the sphere-plate ge- 
ometry can be computed efficiently by Monte Carlo simulations of spin models. Similar to 
the proposal of Hucht fl^] for the plate-plate geometry, we compute differences of the free 
energy by integrating differences of the energy over the inverse temperature. For i? -C Lq, L, 
where R is the radius of the sphere and Lq, L are the linear extensions of the lattice, such 
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TABLE II. Results for the ratio Z^jZ^ of partition functions for the radius R = 3.5. In the 
fourth column we give the number of measurements. In the last column we give the estimate of 
the universal ratio a. In () we give the error due to the statistical error of Z^jZ^ and in [] the 
error due to the uncertainty of Reff- For a discussion see the text. 



h = D + R 


Lq L stat 


Z_ / Zj^ a 


72 


199 160 10^ 


0.1645(28) 


72 


299 250 10^ 


0.1472(26) 


72 


359 300 10^ 


0.1449(26) 


72 


499 400 1.5 X 10^ 


0.1391(21) 7.57(6)[9] 


72 


799 600 10^ 


0.1384(25) 


30 


499 400 1.2 X 10^ 


0.0289(7) 8.43(6) [11] 


40 


499 400 1.2 X 10^ 


0.0552(12) 8.09(6)[10] 


50 


499 400 1.2 X 10^ 


0.0761(15) 8.12(6)[10] 


60 


499 400 10^ 


0.1090(20) 7.71(6) [9] 


90 


499 400 1.2 X 10^ 


0.1830(27) 7.34(6) [9] 



energy differences, when determined naively, are affected by a huge relative variance. Here 
we demonstrate how this problem can be overcome by using the exchange cluster algorithm 



35[. Using the exchange cluster algorithm we define an improved estimator for the energy 
differences with a strongly reduced variance. For a detailed discussion see section [IVl We 
simulated the improved Blume-Capel model, which shares the universality class of the three- 
dimensional Ising model. Improved means that the amplitudes of the leading bulk correction 
are strongly suppressed. We studied strongly symmetry breaking boundary conditions. We 
fixed the spins at the boundary to Sa; = 1. The spins at the surface of the sphere are fixed to 
Sx = s sphere, whcrc we studicd the two cases s sphere = — 1 and 1 in detail. First we verified 
that the method is efficient for a large range of parameters. We studied the whole range 
of temperatures where the thermodynamic Casimir force has a significant amplitude. We 
considered distances between the sphere and the plate up to D^ax = 32, 40, 32 and 10 for 
the radii R = 3.5, 4.5, 7.5 and 15.5, respectively. In order to see a good scaling behaviour of 
the data obtained for different radii, we introduced an effective radius Reff of the spheres 
that we determined by using a particular finite size scaling analysis. For details see the 
appendix |Al We obtain reliable results for the scaling functions -ft'+(x. A) and K-{x, A) for 
the whole range of x and 1 ^ A ^ 12, where x = t[D/$^QY^'^ and A = D/R. Smaller values 
of A would be desirable to make better contact with the Derjaguin approximation and larger 
ones to reach the range of validity of the small sphere expansion. For Ssphere = 1 we find 
that for the values of A studied here, the shape of -ft'+(x, A)A^ viewed as a function of x is 
quite similar to that of limA-s.o Kj^{x, A)A^ obtained by using the Derjaguin approximation. 
For Ssphere = — 1 in the high temperature phase still the same observation holds, however 
in the low temperature phase K_(x, A)A^ is strongly suppressed for A ^ 1 compared with 
A — 0. We attribute this behaviour to the fact that the interface between positive and neg- 
ative magnetisation is closely attached to the sphere, in order to minimize the area of the 
interface. In contrast, for the film geometry, the interface is moving quite freely between the 
boundaries. Finally, by using a different cluster algorithm, which is similar to the boundary 
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flip algorithm [36|, |37[ , we computed the difference of reduced free energies F„ — F+ for large 
distances D at the critical point, where the subscript of F gives the sign of Sgphere- Here we 
made contact with the small sphere expansion. Still for A 30 we see a deviation from the 
leading order of the small sphere expansion at the level of about 7%. 

The present work should be seen as a pilot study. It can be improved and extended in 
many respects. First let us mention that we performed a preliminary study for Ssphere = 0, 
which we do not discuss in the main text. This study shows that for Sgphere = the 
performance of the exchange cluster algorithm is worse than for s sphere = ±1, however still 
physically relevant results can be obtained. Based on this experience we are confident that 
the algorithm discussed here can be applied to a variety of experimentally relevant situations. 
E.g. the crossover of boundary universality classes or patterned surfaces. In particular also 
the sphere-sphere geometry should be accessible by the algorithm discussed here. One might 
also check whether the exchange cluster algorithm allows to improve the results for the film 
geometry that have been obtained so far using standard simulation algorithms. 

One could improve the present study by various means. In order to reach larger lattice 
sizes one might implement the computer program by using the multispin-coding technique. 
To this end one might also try to simulate with a varying lattice resolution: spins far away 
from the sphere could be replaced by block-spins. In order to reduce corrections to scaling 
one might try to improve the spherical shape of the sphere by tuning the couplings at the 
surface of the sphere. 
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Appendix A: Determination of the efTective radius Reff 

In order to determine the effective radius of the sphere, we studied the ratio of partition 
functions r = Z_/Z^ for a lattice of the linear sizes Lq and Li = L2 = L, where the sphere 
is placed exactly in the middle of the system; i.e. h = (Lq + l)/2 for odd values of Lq. The 
spins at Xo = and Xq = Lq + 1 are fixed to s^; = 1. The index of Z refers to the sign of 
Ssphere- At the Critical point, in the scaling limit 

r{R,LQ,L)=G{R/LQ,L/LQ) . (Al) 

In the following, we restrict ourselves to L/Lq = 1/2 and define g{R/ Lq) = G{R/Lq, 1/2). 
Furthermore we require that the ratio of partition functions assumes a fixed value. Our 
ad hoc choice is r = 0.1. For a given radius R this fixes the linear size of the lattice to 
Lq = Lqjix- In order to reduce corrections to scaling we replace Lq by i^o,e// = Lq + Ls 
using Ls = 1.92 jloj. Since L assumes integer and Lq odd integer values, we interpolate 
or extrapolate r{R,LQ,L) to r(i?, Lq, -Lo,e///2). To this end, we simulated for R = 3.5 
and Lq = 131 the linear extensions L = 64, 66 and 68. The results are r = 0.10405(23), 
0.10015(16) and 0.09641(22), respectively. We performed 10'', 2 x 10^ and 10'' measurements 
for L = 64, 66 and 68, respectively. In order to apply the result also to other radii, we 
compute 

- =-0.254(11) (A2) 

0,e//) L/Lo,eff=0.5,Lo=Lo,fi^ 



d{L/L 
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TABLE III. Estimates of r(R, Lq, L^^eff f^) for Lq such that r ~ 0.1. In the last column we give 
the solution of r{R, Lq, LQ^^ff f^) = O-l- Por a discussion see the text. 



3.5 131 0.09927(20) 133.8(3) 
4.5 183 0.09972(30) 185.4(6) 
5.5 235 0.09922(32) 238.6(0.8) 
6.5 287 0.10162(32) 284.7(1.3) 
7.5 331 0.09903(39) 335.9(1.5) 
9.5 435 0.10027(41) 435.8(1.7) 
11.5 535 0.10103(43) 531.9(2.6) 
15.5 731 0.10065(57) 728.6(4.2) 



where we have used the finite difference of the estimates for L = 64 and 68. For larger radii 
we have simulated L = {Lq + l)/2 and have extrapolated r to L = {Lq + 1.92)/2 by using 
eq. (IA2p . In addition to R = 3.5 we have simulated the radii R = 4.5, 5.5, 6.5, 7.5, 9.5, 12.5 
and 15.5. We performed 8 x 10^ measurements for R = 4.5 down to 3.4 x 10^ measurements 
for R = 15.5. In total these simulations took about 5 years of CPU-time on a single core 
of a Quad-Core AMD Opteron(tm) 2378 CPU running at 2.4 GHz. Our results for Lq such 
that r ^ 0.1 are summarized in table UTTl We linearly extrapolate r{R, Lq, Lo^efff^) in Lq. 
To this end, we use 

= 0.11(1) (A3) 

OLq 

which we have computed by using finite differences. Solving r(i?, Lq, -Lo,e///2) = 0.1 with 
respect to Lq we arrive at the estimate of LQ^ffjix given in the last column of table IIIII 
Now we define the effective radius by 

-^e// = C~^LQ^effJix{R) (A4) 

In order to determine the factor c we fitted our data with the Ansatz 

LQ,effjix{R) = c[R + Rs] . (A5) 

where c and Rg are the parameters of the fit and with 

Lo,effjix{R) = c[R + Rs] X (1 + b[R + R,]-^) (A6) 

with the additional free parameter b. The results of our fits are summarized in table IIV[ As 
our final estimate we take c = 49(1) which is consistent with all results given in table HVl 
Finally in table |V] we summarize the results for the effective radii -Re// for the values of R 
that are studied in the remainder of the paper. 

Appendix B: The sphere flip algorithm 

In the Appendix 1 and section IV Fl we have computed the ratio 

r = 1^ (Bl) 
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TABLE IV. Results of fits witli tlie Ansatze IIA5\A6\) . All data with Rmin < R < 15.5 are taken 
into account. In the case of the Ansatz liA6]) we do not report the estimates of b. 





c 


Rs 


xVd.o.f. 


IA5] 4.5 


49.65(22) 


-0.74(2) 


3.11 


|A5] 5.5 


49.08(28) 


-0.65(4) 


1.19 


lA5] 6.5 


49.47(37) 


-0.72(6) 


0.78 


lA6l 3.5 


49.00(31) 


-0.61(5) 


1.97 


lA6l 4.5 


48.57(41) 


-0.52(8) 


2.07 



TABLE V. Results for the effective radius Reff obtained by using eq. (fA4|) and the estimate 
c = 49(1). 



^ Reff 

3.5 2.73(6) 

4.5 3.78(9) 

7.5 6.86(17) 

15.5 14.87(38) 



by using a special cluster algorithm. The index of the partition function Z refers to the sign 
Ssphere of the Sphere. This cluster algorithm is closely related to the boundary flip algorithm 
discussed in refs. [sgI, 37| that allows to compute the ratio Za/Zp, where p refers to periodic 



and a to anti-periodic boundary conditions. In ref. |23[ we had employed a similar algorithm 
to compute where the indices of Z refer to the signs of the spins at the boundaries. 

Let us consider the composed system 

Z = Z_ + Z+= Yl Yl ^M-Hi{s}, Ssphere)) (B2) 

^sphere {"^l 

where Sgphere can be viewed as an Ising spin. This system can be updated by using the 
standard cluster algorithm. When the cluster that contains the sphere is not frozen to the 
boundary, s sphere can be flipped. In this composed system, the ratio r = is given by 

the ratio of expectation values r = ^ . ) / i^i s ^ )■ 

-C^ \ ^ 'i'^ sphere I I \ ^•i'^ sphere I 

In fact it is sufficient to simulate only. Since for s sphere = ~1 the cluster that contains 
the sphere never freezes to the boundary it is sufficient to compute the fraction of cluster 
updates that would allow to update Ssphere = 1 to Ssphere = —1- To this end, on has to 
check for each cluster update, whether the cluster that contains the sphere is frozen to the 
boundary. In the composed system 

z+ p{ — y +) 

where p{ )■ +) and p( )■ +) are the probabilities to update Ssphere = —1 to Ssphere = 1 

and vice versa, respectively. As discussed above p{ +) = 1. 
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We have updated the Sgphere = 1 system by a combination of the local heat-bath algorithm 
and the Swendsen-Wang cluster algorithm 4^. In the case of the Swendsen-Wang cluster 



algorithm [42] the boundary conditions have to be taken into account: All spins that belong 
to clusters that are frozen to the boundary or contain the sphere keep their value. All other 
clusters are flipped with the probability 1/2, where flipped means that the sign of all spins 
in a given cluster are changed. 
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